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1.0 SUMMARY 

The analysis, structure, and capability of a generalized 
precision interplanetary trajectory computation program are 
discussed, with emphasis being placed on the description of 
input and output. Sample cases showing input and output in- 
formation are included. 

2.0 INTRODUCTION 

This report describes a computer program which was de- 
veloped by the Texas Center for Research of Austin, Texas, for 
computing precision interplanetary trajectories. The program 
uses numerical optimization techniques and a variation of 
Newton's method to produce trajectory iterates which converge 
to an approximate n-body interplanetary trajectory. This ap- 
proximate trajectory is then used with a perturbation procedure 
to produce the precision n-body interplanetary trajectory. 

This procedure was first tested on earth-moon trajectories 
where it was used for a published paper (ref. 1) and several 
reports (refs. 2, 3). 

The procedure developed in this study has been used to 
compute two-dimensional and three-dimensional earth-moon tra- 
jectories and two- and three-dimensional interplanenery n-body 
trajectories using circular planetary orbits, elliptic plane- 
tary orbits, and planetary orbits specified by the JPL Analy- 
tical Ephemeris. In addition, trajectories from earth to Venus 
to a heliocentric radius vector with a Av at Venus have been 
computed. The program was then used in an iterative mode to 
minimize the Av at Venus in order to simulate the free-fall 
fly-by. The modifications necessary to provide this capability 
will be made in the program and an addendum to the User's 
Manual will be provided. 
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This program, possibly with some minor modifications, will 
fill a variety of needs for precision interplanetary trajec- 
tories. This report will not cover all possible modifications 
of the program, but will indicate the types of modifications 
which could be made and will outline the potential uses of 
various modified versions of the program. 


3.0 PROBLEM ANALYSIS 

The trajectory computation problem for free-fall trajec- 
tories can be stated as follows: 

Find the set of initial velocities v(t^) such that a 
particle subject to the equations of motion 


moves from a known initial point (x^,t^) to a known final 
point (x^jt^) • This is a classical two-point boundary value 
problem. It is not an optimal control problem as there are no 
controls in this formulation. 

The details of the derivation will not be covered here. 
Only the important equations and how the technique works will 
be presented. 


Primary Convergence Algorithm 


Rewrite the equations as 


x = u(t) 

A = F (x , t) 


C 2 ) 


with initial conditions A ( t ^ ) = u(t^) . Clearly, we want to 
change the system in such a way that u(t) and A (t) approach 
the same values, i.e., u(t) . Note that if a u(t) history 
is guessed (a control program in the optimal control vocabulary) 
then a position history, x(t) , and a history of the variables 
A(t) can be obtained. 
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The procedure used to drive u(t) and A(t) to the same 
values requires that both of equations (2) above be linearized. 
The new values for u(t) and A(t) are given by 


u new (t) = u(t) + 6u(t) 
A new (t) = + 6A(t) 


(3) 


Setting u new ( t ) equal to > one obtains 

Su(t) = 6A(t) + A(t) - -u(t) 


(4) 


where u(t) and A(t) are evaluated on the previous iterate. 

In order to control stepsize and to keep consecutive iter- 
ates close enough to each other so that linearization is valid, 
the scalar variable P(0 <_ P _< 1) is introduced into equation 
(4) . Equation (4) is modified in the manner shown below 

fiu(t) = 6 A (t) + P [A (t) - u(t)] (5) 


Linearization of equations (2) leads to 

Sx = u(t) 

6A = F x 6x(t) 


( 6 ) 


where is a 3 x 3 matrix of partial derivatives of the 

F 's with respect to the position coordinates. Substitution 
for 6u in equation (6) from equation (5) leads to 

Sx = 6A(t) + P[A(t) - u(t)] 

SA = F 6x(t) 
x K J 

It is now assumed that both <5x(t) and SA(t) are linear 
functions of the initial values of Sx^ and SA^ . Since 
J x. = 0 , it is possible to write 
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6x(t) = A(t) 6A ± + M (t) 
6A(t) = B ( t) 8X i + N (t ) 


( 8 ) 


where A and B are 3x3 time -dependent matrices and M(t) 
and N(t) are time- dependent 3-vectors. The conditions which 
define A, B, M and N will be determined below. 

From the conditions ox(t^) = 0 and 6A(t^) = 6A^ , the 
initial conditions on A, B, M and N can be defined; i.e., 

A(t i ) = 0 

B(t.) = I (the 3x3 identity matrix) 

1 (?) 

M(t i ) = 0 

N(t i ) = 0 

Taking derivatives of equations (8) and substitution of 
these results plus those of equations (8) into equations (7) 
leads to the relations 


A6A i + M = B6A i + N + P(A - u) 
B <5A ■ + N = F ASA. + F M 

1 _ X X X 


( 10 ) 


Thus, the conditions defining A, B, M and N are 


A = B 

A(t t ) 

= 0 


- B = F x A 

B(t i ) 

= I 

M = N + P (A - u) 

M(t i ) 

= 0 

• 

N = F M 

X 

N(t i ) 

= 0 


From these relations, A(t), B(t), M(t) , and N(t) can be 
determined. Then, 

SA i = A~ 1 (t f ) [ 6x (t £ ) - M(t f ) ] 


( 12 ) 
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and 


6A(t) = B(t)A~ 1 (t £ ) [6x(t £ ) - M(t f )] + N (t ) 


(13) 


Also , 


6u(t) = B(t)A" 1 (t f ) [(5x(t £ ) - M(t £ ) ] + N(t) + P[A(t) - u(t) ] (14) 


The computational procedure is as follows: 

1. Guess u(t) 

2. Integrate the equations below from t- to t £ with 


the indicated initial conditions: 
x = u ( t ) 

A = F(x,t) 

A = B 


x(t i ) = x i , given 
A(t £ ) = u(t i ) 


B = F x A 


M = N + P(A - u) 


N = F M 

X 


A(t.) = 0 
B( ti ) - I 
M(t.) = 0 
N (t £ ) = 0 


(15) 


3. Using the specified values for x £ and the integrated 
values of x(t £ ) , calculate 


6x £ = x £ - x ( t £ ) 


(16) 


4. Next, calculate 


6A i = A' (t £ ) [ 6x f - M f ] 


(17) 


and 


- i 


6u(t) = B (t) A • (t f )[6x f - M f ] + N(t) + P(A(t) - u(t)] 


(18) 


5. Now set 


A . 


= A . + 6A . 


-new 


u new (t) = + 


( 19 ) 
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6. Return to step 2 and repeat the process unti l A(t) 
is sufficiently close to u(t) . 

Once the process described in steps 1 through 6 above 
converges, the values (x(t^),A^) are used as initial condi- 
tions (v(t^) = A^) in the integration 

x = v x. SDecified 

l 

v = F (x , t) v ( t i ) = X ± 

This integration is employed in a standard perturbation proce- 
dure to refine the initial values v(t^) (i.e., to refine A^ ). 

The perturbation procedure normally converges in three or four 
iterates. Note that the standard perturbation algorithm is em- 
bedded in the computational algorithm being used. As A(t)-> u(t) , 
the term P(A - u) goes to zero and M and N become identi- 
cally zero. When this occurs, equation (12) becomes the stan- 
dard equation for SA- from the perturbation algorithm and 
equations (13) and (14) become the same. 

4.0 PROGRAM STRUCTURE 

The n-body interplanetary trajectory program uses four 
basic groups of subroutines. The functions of the four groups 
of subroutines, and the individual subroutine names are listed 
in Table I . 

The MAIN Program calls RDATA for input (see Table II for 
input structure) and then converts the planetocentric equa- 
torial inputs to heliocentric ecliptic coordinates using sub- 
routines COORDS and EPHEMP . Next MAIN calls subroutine DRIVER 
which controls the program from this point forward. 

Subroutine DRIVER sets up the initial conditions for the 
first iteration and either generates a u(t) history from a 
heliocentric trajectory patched to hyperbolic orbits near the 
planets (subroutine VGUESS) or generates a u(t) history by 
integrating a guessed set of initial velocities for the trans- 
fer in the n-body solar system employing RK 3(4). 
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Once the initial u(t) history is known, DRIVER controls 
the iteration described under PROGRAM ANALYSIS until u(t) 
and A ( t ) match sufficiently well to allow the shift to a 
standard perturbation procedure for final convergence. This 
iteration process employs RK 3(4) for integration and SPLINE 
for interpolation. 

When the shift is made to the perturbation procedure, the 
RK 7(8) integrator is employed. Both RK 3(4) and RK 7(8) 
obtain derivatives from subroutine DERIV. 

Slightly different derivatives are called for by each 
procedure, so DERIV contains a branch. Subroutine DERIV uses 
subroutine EPHEMP and SPLINE to aid in the generation of the 
required derivatives. 

5.0 PROGRAM CAPABILITY 

The program possesses four options which may be used as 
described below. The option chosen is indicated by the input 
parameter N0PT (N0PT=1,2,3, or 4). The options and their uses 
are described below. 

Option 1 (N0PT=1) 

When operating in this mode, the program finds the n-body 
trajectory which goes from (t^,x^) to (t^,x^) where x^ is 
the specified initial position vector, t^ is the initial time, 
Xf is the specified final position vector, and t^ is the 
final time. 

The initial and final position vectors a r e normally input 

in kilometers in planetocentric equatorial coordinates with x^ 
being measured in a system centered at planet NP(1) and X£ 
measured in a system centered at planet NP(2). The planetary 
indices used are Mercury=l, Venus=2, Earth=3, etc. If it is 
desired to input either x^ or x^ in heliocentric ecliptic 
coordinates in units of AU, the corresponding planetary index, 
NP(1) or NP(2), should be set to 10. 

For this option, the program converts the initial and 
final position vectors into heliocentric ecliptic coordinates 
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(if necessary) , 

and then 

obtains an initial velocity history 

from a Lambert 

solution 

to start the first order iteration 

process . 



For N0PT=1 

runs, the required input data is: 

TP (1) 

t i ’ 

the initial time (J.D.) 

TP ( 2) 

t f ’ 

the final time (J.D.) 

TP (3) 


unspecified 

XI(1) 
XI (2) 
XI (3) 

x i ’ 

the initial position vector in km. if 
NP (1 ) <9 in AU if NP(1) > 10 

XF(1) 
XF (2) 
XF ( 3) 

x f , 

the final position vector in km if 
NP (1 ) < 9 in AU if NP(1) > 10 

VI(1) 


unspecified 

VI (2) 



VI (3) 



VF (1) 


unspecified 

VF (2) 



VF (3) 



NP ( 1 ) 


initial planet index 

NP (2) 


final planet index 

N0PT 


integer, =1 

NPRT 


integer, print key 


NPT=0, do not print last iterate 
NPRT=1, print every NSPRT^h integra- 
tion steps on last iterate 
(converged iterate). 

NP_RT_= t_l, p.ri n t_ every- . NSRRTib 

tion on every iterate. 

NSPRT integer, see NPRT 

Options' 2 and 3 (N0PT=2/N0PT=3) 

These options are available in order to allow a run which 
has been terminated prematurely (due to time limit, for example) 
to be restarted at the last iterate computed. If the program 
was still operating in the first order iteration mode when the 
termination occurred, N0PT=2 should be used. If, however, the 
shift to the perturbation method has already occurred, N0PT=3 
should be used. 
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In cither case, the initial position vector x , (XI) and 

the initial velocity vector v- (VI) will be in heliocentric 


eel i pt ic 

coordinates in 

units of AU and AU/year. These values 

should be 

: used directly 

as input for the N0PT=2/N0PT=3 options. 

For these 

; options 

t- 

1 

TP (1) and t £ TP(2) as in N0PT=1. 

For 

N0PT=2/N0PT=3 

runs, the required input data is: 


TP (1) 

*i > 

initial time 


TP (2) 

t £ , 

final time 


TP (3) 


unspecified 


XI (1) 


initial position vector in AU 


XI (2) 

X i > 



XI (3) 




XF (1) 


final position vector in AU 


XF (2) 

x £ , 



XF (3) 




VI (1) 


initial velocity vector in AU/yr. 


VI (2) 

V i ’ 



VI (3) 




VF (1 ) 


unspecified 


VF ( 2) 




VF (3) 




NP (1) 


unspecified (can be specified but will 


NP (2) 


not be used) 


N0PT 


integer, =2 or 3 


NPRT 


as in N0PT=1 option 


NSPRT 



Option 4 

(N0PT=4) 




The purpose of this option is to enable the user to compute 

a trajectory which goes from (x^, t £ ) through a point (x^, t m ) 

to a point (x f , t y ) with a velocity impulse, AV , at (x m , t m ) 

and then, by allowing the point x m to vary, make the AV at 

(x , t ) as small as possible. This option uses information 
v m’ m' 

obtained from two N0PT=1 runs as input data. 

In order to obtain the information required for input for 

the N0PT=4 option, the user must set up an N0PT=1 run from 

(x- , t.) to (x , t ) and get the converged initial position 
v i ’ x v m m 
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and velocity components at t^ (in heliocentric ecliptic co- 
ordinates) . Then he must set up an N0PT=1 run from (x , t ) 

m m' 

to (*£, t £ ) and get the converged final position and velocity 
components at t£ (in heliocentric ecliptic coordinates) . 

These values are then used as XI, VI, XF , and VF (all 3-vectors) 
in the input for the N0PT=4 run. The times are input in the 
following manner: t^ TP(1), t TP(2), t ( TP(3). The 

X la 1 

quantities NP(1) and NP(2) do not need to be specified for an 
N0PT=4 run. 

For N0PT=4 runs, the required input is: 


TP (1) 

t i ’ 

initial time (J.D.) 


TP (2) 

t m ’ 

fly-by time (J.D.) 


TP (3) 

t f ’ 

final time (J.D.) 


XI (1) 
XI (2) 

x i > 

initial position vector 
N0PT=1 run) 

in AU (from 

XI (3) 
XF (1 ) 
XF (2) 

x f , 

final position vector in 
N0PT=1 run) 

AU (from 

XF (3) 
VI (1) 
VI (2) 
VI (3) 

v i ’ 

initial velocity vector 
N0PT=1 run) 

in AU/yr. (from 

VF (1) 
VF (2) . 
VF (3) 

v f ’ 

final velocity vector in 
N0PT=1 run) 

AU/yr. (from 

-N-P-C-1-) 


unspecified 


NP (2) 




N0PT 


integer, =4 


NPRT 


as in N0PT=1 


NSPRT 






TABLE I 


Basic Subroutine Groupings 
for the 

Interplanetary Trajectory Generation Program 

1. Control Routines 

A A A T At 
! v ii V. X IN 

DRIVER 

2. Numerical Integration Related Routines 

RK34C0 

RK34 

RK78C0 

RK7 8 

DERIV 

VGUESS 

SPLINE 

3. Celestial Mechanics Routines 
C0NIC 
0ET0RC 
EPHEMP 
TRFM 
C00RDS 
EULER 
DLMBRT 
UVFNS 

4. Vector Manipulation Routines 

UNIT 

D0T 

VMAG 

CR0SS 
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TABLE II 


Description of Input 


Variables 

Description 

T P ( 1 ) , T P ( 2 ) , TP (3) 
Units (Julian days) 
Double Precision 

Times in Julian days as required for the 
various problem options 

XI (1) ,XI (2) ,XI (3) 
Units (km or AU) 
Double Precision 

Initial position vector for N0PT=1, XI is 
expressed in planetocentric equatorial co- 
ordinates at planet NP(1) in kilometers for 
N0PT=2,3, or 4, XI is expressed in helio- 
centric ecliptic coordinates in A.U. 

XF (1) ,XF(2) , XF (3) 
Units (km or AU) 
Double Precision 

Final position vector (see XI for details 
on units for various values of N0PT) 

VI (1) ,VI(2) ,VI (3) 
Units (AU/yr) 
Double Precision 

Initial velocity vector in heliocentric 
ecliptic coordinates in AU/yr. (Needed 
only for N0PT=2,3, or 4) 

VF (1) , VF ( 2 ) , VF (3) 
Units (AU/yr) 
Double Precision 

Final velocity vector in heliocentric 
ecliptic coordinates in AU/yr. (Needed 
only for N0PT=4) 

NP ( 1) ,NP (2) 
Integers 

Planet index denoting departure planet and 
target planet in system where Mercury=l, 
Venus=2, etc. 

N0PT 

Integer 

Option choice key, N0PT=1,2,3, or 4. 

NPRT 

Integer 

Print option key 

NPRT=-1 Print every NSPRT r integration 


step on every iterate 
NPRT=0 Print only summary information- - 
do not print at intervals on any 
iterate , 

NPRT=1 Print every NSPRT r integration 

~~ — — — — s-t-e.p_.Q-n_t.h_e_. l ast f conve rged ) 

iterate only 

NSPRT The number of integration steps between 

Integer prints if NPRT/O. (see NPRT) 


For more detail on input, see Section 5.0, PROGRAM CAPA- 
BILITY. 

The data is input using the system routine RDATA. The call 
to RDATA is 

CALL RDATA (TP , XI ,XF,VI , VF , NP ,N0PT , NPRT , NSPRT) 
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Output Quantity 

ITERATION 

(Integer) 

FINAL TIME 
Units (years) 

INTEGRATION STEPS 
(Integer) 

XI 

Units (AU) 

XF 

Units (AU) 

VI 

Units (AU/yr) 

VF 

Units (AU/yr) 

DELXF 
Units (AU) 


DUMAX 

Units (AU/yr) 


DELVI 

Units (AU/yr) 

MINIMUM 

APPROACH 

I) I STANCES 
Units (AU) 

F [NAL 
APPROACH 
DISTANCES 
Units (AU) 


TAP. IP. I I I 

Description o P Output 


Description 

A counter which shows the current number 
of iterations used during any one mode 
of operation 

Elapsed time during the current iteration 
from the initial time 

Total number of steps taken by the inte- 
grator to complete the current iteration 

Initial position in heliocentric ecliptic 
coordinates (3 components) 

Final integrated position in heliocentric 
ecliptic coordinates (3 components) 

Initial velocity (3 components) 

Final integrated velocity (3 components) 

Terminal position miss defined as differ- 
ence between input final position and 
integrated final position (3 components) 

Maximum change in the control velocity 
history for the next iteration of this 
mode 1 operation. This quantity is mean- 
ingless during modes 2 and 3 operation 
(3 components) 

Change in the initial velocity for the 
next iteration (3 components) 

The minimum radial distance obtained any- 
ti me durin g the cur rent iteration from 
each of the solar system memb e r s FisTecl ' 

The final radial distance obtained during 
the current iteration from each of the 
solar system members listed 



fraction of planets g « a v t T r used 
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TABLE V 

Sample Data Set for N0PT=4 Case 

The following data set will be assumed to be punched on 
cards started in column 1 with no spaces unless indicated. 

Card 1--TP 

Dl=2. 4443299842 D+ 06 , 2 . 4445049474D+06 , 2 . 4446221842D+06 
Card 2 - - XI 

D2= - 9 . 824794 567024D-01 , -1 . 8206340107 50D- 01 , 9 . 417 4 1 311 3088D- 06 
Card 3--XF 

D3=2.610368798355D-01 , - 4 . 199 89 7439 846D- 01 , - 3 . 679 84643701 0D- 02 
Card 4- -VI 

D4=2. 691275379500 D+ 00, -4. 440086 5 78800D+ 00,-1. 32295707370 0D -01 
Card 5--VF 

D5=8. 7038421339 D+ 00, 4. 30945631490D+00, -6. 446741114700D-01 

Card 6--NP 
16=10,10 

Card 7--N0PT 

17 = 4 

Card 8--NPRT 

18 = 1 

Card 9--NSPRT 

19 = 10$ 

Note that the dollar sign ($) should follow immediately the 
last digit of data in each data set in RDATA input. 



